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Abstract 

We consider the problem of reconstructing general solutions to the Helmholtz equation Au + X 2 u = 0, 
for some fixed A > 0, on some domain U C R 2 from the data of these functions at scattered points 
xi, . . . ,x n C ft. This problem typically arises when sampling acoustic fields with n microphones for the 
purpose of reconstructing this field over a region of interest ft that is contained in a larger domain D over 
which the acoustic field is generated. In many applied settings, the boundary conditions satisfied by the 
acoustic field on dD are unknown as well as the exact shape of D. Our reconstruction method is based 
on the approximation of a general solution u by linear combinations of Fourier-Bessel functions or plane 
waves ek(x) := e lk a: with jkj = A. We study two different ways of discretizing the infinite dimensional 
space V\ of solutions to the Helmholtz equation, leading to two different types of finite dimensional 
approximation subspaces, and we analyze the convergence of the least squares estimates to u in these 
subspaces based on the samples (u(xi))i=i,..., n - Our analysis describes the amount of regularization that 
is needed to guarantee the convergence of the least squares estimate towards it, in terms of a condition 
that depends on the dimension of the approximation subspace and the sample size n. This condition also 
involves the distribution of the samples and reveals the advantage of using non-uniform distributions that 
have more points near or on the boundary of fi. Numerical illustrations show that our approach compares 
favorably with reconstruction methods using other basis functions, and other types of regularization. 

1 Introduction 

A common problem in acoustics is to obtain a precise approximation of the soundfield over a spatial domain 
ft of interest, using the smallest possible number of pointwise measurements, e.g. as provided by micro- 
phones. For instance, one may wish to measure the complex radiation pattern of an extended source (source 
identification problem), to localize a number of point sources within a spatial domain (source localization 
problem), or to optimize the output of a sound reproduction system over a large control area, to name only 
a few applications. In practice, the main difficulty that one is usually faced with is how to handle reverber- 
ation : the reverberant field might well be of a magnitude comparable to the direct sound, and it depends 
in a non-trivial way on both the geometry of the domain D where the acoustic field is defined and the type 
of boundary conditions on dD (with Dirichlet or Neumann as ideal cases, but more likely in engineering 
problems with a frequency-dependent mixed behaviour). 

The goal of this paper is to study the accuracy that can be achieved when approximating the acoustic field 
over the domain ft C D, based on a set of point measurements, without precise knowledge on the geometry 
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of D and boundary conditions on dD. A general setting is the following: the soundfield p(x,t) is measured 
at microphones located at positions Xi, . . . ,x n € f2, and over a (discretized) time interval [0, T]. After 
application of the (discrete) Fourier transform J- in the time variable, and considering a given frequency uj, 
the function 

u{x) := J-p(x, u) 

is a solution on D to the Helmholtz equation 

Am + A 2 u = 0, (1.1) 

where A = ui/c, with c denoting the wave velocity, and where the boundary conditions are unknown to us. 

Depending on the applications, the geometry of the domain f2 may either be 2-D (membranes) or 3-D 
(rooms). Here, we restrict our attention to the 2-D setting, although certain features of our approach are 
readily extendable to the 3-D setting. Our problem therefore amounts to reconstructing, on some domain 



£1 C R , a general solution to the Helmholtz equation (1.1 1 from its sampling at points Xi,...,x n € fl 



These samples may be measured exactly or up to some additive noise. We denote by 

yi = u{x{) + -qi, l = l,...,n, (1.2) 

these samples, where rji represent the additive noise. 

Reconstruction from scattered points is a widely studied topic, and a variety of methods have been 
proposed and analyzed. Many existing methods can be viewed as reconstructing some form of approximation 
to the unknown function u by simpler functions such as splines, partial Fourier sums or radial basis functions. 
The success of these methods therefore relies in good part on the quality of the approximation of u by such 
simpler functions, which is typically governed by the smoothness of u. 

In our present setting, the fact that u obeys the Helmholtz equation, may be used in addition to its 
smoothness in order to guarantee the accuracy of certain approximation schemes, which are well adapted to 
such solutions. We introduce the Fourier-Bessel functions 

& AJ (z):=e^J 3 -(Ar) (1.3) 

where (r, 9) are the polar coordinates of x and Jj is the j-th Bessel function of the first kind. b\ j is solution 



to the Helmholtz equation ( 1.1 1 over K. if and only if its parameter A is the same as in ( 1.1 1. Denoting V\ 



the set of the solutions of ( 1.1 1, it is known jjj that 



-L 2 (0) 



V x = span{6 AJ } (1.4) 



and that the solutions of (1.1) can be approximated by elements of the subspaces = span{6 A ,ii — m < 
j < m} as m grows. 

An alternative approximation scheme uses planes waves defined by 

e k (x) := e ik -* (1.5) 

which are solutions of ( |1.1| if and only if |k| = A. The spaces V^j, spanned by the particular plane waves 



ekj ! kj 



can also be used to approximate solutions of (1.1 ) as m grows [3 
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The most widely used approach to approximate u in a finite dimensional space V m , from its data at 
points Xi, ...,x n , is the least squares method, namely with m < n solving the minimization problem 



l j2\ yi -v{ Xl )\ 2 . (1.7) 



7r = argmin„ eK 

»=i 

The effectiveness of the least squares approximation is governed by a certain trade-off in the choice of the 
dimension m of the approximation: 

• A small value of m leads to a highly regularized reconstruction of u, which is usually robust but has 
poor accuracy. 

• A large value of m may lead to unstable and therefore inaccurate reconstructions although the space 
V m contains finer approximants to u. 

Let us observe that regularization is relevant even in a noiseless context where the function is measured 
exactly: for example choosing m — n corresponds to searching for an exact interpolation of the data which 
may be very unstable and inaccurate, similar to the Runge phenomenon in polynomial approximation. 

In this paper, we discuss the amount of regularization which is needed when applying the least squares 
method using the finite dimensional subspaces and extracted from V\. With such discretizations, the 
distribution of the sampling points x%,...,x n has an influence on the above described trade-off. Our main 
theoretical result, established in the case of a disc, shows that higher values of m, leading therefore to better 
accuracy, can be used if the Xi are not uniformly distributed on f2 in the sense that a fixed fraction of these 
points are located on the boundary dfl. This result is confirmed by numerical experiments. 

The rest of this paper is organized as follow: we give a brief account in §2 on approximation of solutions to 



1 1.1 1 by Fourier-Bessel functions and plane waves which relies on Vekua's theory, and in §3 on general results 
on the stability and accuracy of least-squares approximations recently established in [2] . We then study in §4 
the spaces and Vj^ in more detail, in the particular case where f2 is a disc, and use the above mentioned 
results to compare least-squares approximations on these spaces based on different sampling strategies. In 
§5, we present numerical tests that illustrate the validity of this comparison. We also show that our approach 
compares favorably with reconstructions based on other approximation schemes such as partial Fourier sums 



(that do not exploit the fact that it is a solution to (1.1)) and to other form of regularizations such as 
Orthogonal Matching Pursuit. 



2 Approximation by Fourier-Bessel functions and plane waves 



Results on the approximation of solutions to (1.1) by Fourier-Bessel functions and plane waves given in 
[4] are based on the theory developed in the 1950's by Vekua [6 . This theory generalizes approximation 
results for holomorphic functions, viewed as solutions of Au = 0, to solutions of more general elliptic partial 
differential equations, by means of appropriate operators that link the two types of solutions. 

In the case of the Helmholtz equation on a domain ft, that is star-shaped with respect to a point which 



is fixed as the origin 0, these operators (in their version mapping harmonic functions to solutions of (1.1)) 
have the explicit expression 

X\x\ f 1 



Vi<f>(x) = --^-J : ^=Jx{\\x\Vl^t)dt, (2.1) 
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and 

1 

V 2 cf>(x) = J>{x) + MJ ^_}_^ h (X\x\VT^t)dt, (2.2) 

where \x\ stands for the euclidean norm of ir, J\ is the Bessel function of the first kind of order 1 and I\ the 
modified Bessel function of the first kind of order 1, see [3] for more details. These operators have important 
properties: 

• Their are linear. 

• Vi maps harmonic functions to solutions of the Helmholtz equation, and V2 does the converse. 

• They are continuous in the Sobolev H k norms for all k > 0. 

• They are inverse to each other on these spaces. 

As a consequence, any approximation method for harmonic functions can be translated as an approxi- 
mation method for solutions of the Helmholtz equation. In particular, approximation of harmonic functions 
by harmonic polynomials of degree m translates as approximation of solution of the Helmholtz equation by 
the so-called generalized harmonic polynomials m which are their image by V±. The generalized harmonic 



polynomials of degree m can be expressed as linear combinations of the Fourier-Bcssel functions (1.3|, leading 
therefore to results for the approximation of solutions to (1.1) by elements of in Sobolev norms. More 



precisely, the following result can be obtained when the domain Q is convex, see theorem 3.2 of [1] 

log m x p 



min \\u — v\\ H k < C IMIi?"; (2-3) 

i>ev£ \ m J 

where the constant C depends on p, k, A and the geometry of Q. This results still holds for more general, 
star-shaped convex domains, with a slower convergence. 

Planes waves and Bessel functions are related by the Jacobi- Anger identity 

e k, = ]Tz m J m (Ar)e m ^. (2.4) 
where := A(cos(</>), sin(0)), and its converse, the Bessel integral 

7T 

U\r)e me = / e^e in <t>dd>, (2.5) 



Approximating the integral in ( 2.5 ) by a discrete sum, with uniformly sampling of the wave vectors k^ on the 



circle of diameter A, leads to approximations of solutions to ( 1.1 ) by linear combinations of the 2m + 1 planes 
waves ekj, that is, by elements of Vf n . It is also known (see theorem 5.2 of [1]) that such approximations 
have the same convergence properties, e.g., for a convex domain, 

min ||«-«|| H * < C { J \\u\\ HP . (2.6) 

vevi> \ m I 
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3 Least-square approximations 



The results of the previous section quantify how a general solution u of the Hclmholtz equation can be 
approximated by functions from spaces or V^. We are now interested in understanding the quality of 
approximations from these spaces built by the least squares methods based on scattered data (xi, yi)i=i,...,n- 
In particular, we want to understand the trade-off between the dimension m and the number of samples 



n. Ideally we would like to choose m large in order to benefit of the approximation properties (2.3 1 and 



(2.6), however not too large so that stability of the least-square method is ensured. We are also interested 
in understanding how the spatial distribution of the sample x\ influences this trade-off. 

This problem was recently studied in [2J, in a general setting where the xi are independently drawn ac- 
cording to a given probability measure v defined on Q. This measure therefore reflects the spatial distribution 
of the samples. For example, the uniform measure 

dv := \Q\- X dx, (3.1) 

tends to generate uniformly spaced samples. In order to present the general result of [2], we assume that 
(V m ) m >i is an arbitrary sequence of finite dimensional spaces of functions defined on Q with dim(y m ) = m. 
We introduce the L 2 norm with respect to the measure v 

1/2 



|M|:=(/mW) , (3.2) 



and we define the best approximation error for a function u in this norm as 

a m (u) := min \\u — v\\ (3-3) 

Note that, in the noiseless case, the least squares method amounts to computing the best approximation of 
u onto V m with respect to the norm 

ihi« : =(^Ek^)I 2 ) 1/2 - ( 3 - 4 ) 

This norm can be viewed as an approximation of the norm \\v\\ based on the draw, and it is therefore natural 



to compare the error ||u — tt|| where 7r is computed by (1.7) with a m (u). We give below a criterion that 
describes under which condition on m these two quantities are of comparable size. 

Here, we assume that (Li, . . . ,L m ) is a basis of V m which is orthonormal in L 2 (fl,i/). We define the 
quantity 

m 

K{m) = K(V m , v) := maxV |L,(x)| 2 , (3.5) 

i=i 

which depends both on V m and on the chosen measure v, but not on the choice of the orthonormal basis 
since it is invariant by rotation. 

We also assume that an a-priori bound ||u||^oo < M is known on the function u. We can therefore only 
improve the least squares estimate by defining 

u:=T M (ir), (3.6) 



where 7m(£) := sign(t) max{|i|, M} and ir is given by (1.7). The following result was established in [2], in 



the case of noiseless data, i.e. r\i — in (1.2) 
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Theorem 3.1 Let r > be arbitrary but fixed and let n := ^tt§~ ■ If m is such that 

n 

K(m)<Kz , (3.7) 

log n 

then, one has 

E{\\u - u|| 2 ) < (1 + e{n))a m {u) 2 + 8M 2 n~ r , (3.8) 
where e(n) := — > asn-> +00. 



It is also established in [2] that the condition (3.7) ensures the numerical stability of the least-square 
method, with probability larger than 1 — 2n~ r . These results suggest to set the regularization level by picking 
the largest value of m* — m*(n) such that (3.7 1 holds. The dependance of m*(n) with n is obviously related 
to that of K(m) with m. In particular, slower growth of K(m) with m implies faster growth of m*(n) and 
therefore faster convergence of the least squares approximation. Notice that we always have 

/m 
\ L i? dv = m - (3-9) 

In the next section, we evaluate K(V^, v) and K(V£, v) in the case where is a disk, for various choices of 
the measure v. 



4 Fourier-Bessel functions and plane waves on the disc 

As mentioned in the previous section, the quantity K(m) depends both on the space V m and the measure 
v that reflects the sampling strategy. Here we study the case where V m is either one of the spaces of plane 
waves or of Fourier-Bessel functions defined in the introduction. Note that the dimension of these 
spaces is not m but 2m + 1, but for notational simplicity we maintain the shorthand K(m) for K(V^,u) 
and K{V§,v). 

We consider two sampling strategies. The first one uses the uniform probability distribution 

Vo = WY ( } 

and the second one combines uniform sampling on the domain and on its boundary, with proportion < 
a < 1, according to the probability distribution 

% dx da , , „. 

^ : - {1 ~ a) w\ +a w\- (4 - 2) 

In order to obtain explicit results, we focus on the simple case where 17 is a disk. Without loss of generality, 
we fix 

:= {x e W 2 : \x\ < 1}. (4.3) 
4.1 Fourier-Bessel approximation 

Fourier-Bessel functions are orthogonal on the disk with respect to any rotation-invariant measure, because 
of their angular dependence in e m9 . This allows a simple computation of the quantity K(m) for the space 
V^, leading to the following result. 



G 



Theorem 4.1 For the space V^ n and the measure v a on the unit disk ft, one has for sufficiently large m 

K(m) > c + cirn 2 , (4.4) 

when a = (that is, for the uniform measure), for any c\ < 1/4, and where cq depends on c\ and X, and 

2m 



K(m) < C + 

a 

when a > 0, where C depends on A and a. 

Proof: Since the Fourier-Bessel functions are orthogonal in L 2 (Vl,v a ), we have 



K(m) 



E 

j=-m 



IM 2 



L~(fi) 



In the case of the uniform measure, we bound K(m) from below. We first write 



K(m) > 



E 

j=-m 



L~(an) 



We next bound \\bj\\ 2 , for j > [A] (the case j < — |~A] is identical, as \bj \ = \b-j\), according to 



^-(Alxl) 2 ^ 



n 
i 



2 / rJj{\r) 2 dr 



(4.5) 



(4.6) 



(4.7) 



A 2 



^0' + l + 2p)J? +1+2p (A) 



p=0 
oo 



^ ^E(j + 1 + 2 p)( a /j) 2+4p ^(a) 

,3 + 1 ( 



1 



(A/j) 4 



f \i-(Vj) 4 .7 + 1 (1 - (A/j) 4 ) 



2 J 



J 7 2 (A) 



where the third equality is identity (11.3.32) of pQ for the third equality, and the first inequality comes from 
(A.6) of 0. As |6j(l,6>)| = |Jj(A)|, we have 



> 



j 



i 



(A/i) 4 



INI 2 ~ 4(j + 1) \1 - (A/j) 4 J + 1 (1 - (A/j)4) 2 / 
For any c < 1, there exists a Jo such that when j > jo, 

MM) 2 U J 



Nl 2 " c 4 



and 



3 = -Jo 



l^(A)| ; 



2 ^ E 



j <i<m 



(4.8) 

(4.9) 

(4.10) 
(4.11) 



7 



which proves the bound (4.4 1 when m > jo 



In the case of mixed sampling, we bound K(m) from above by 

K(m) < f] lMl " {Q) 



IM 2 



When j > A, the function r i-> J/(Ar) is monotone increasing on [0, 1], so that = Jj(^)- Thus, 



0jll£°°(fi) 



2jt 



r> o 



< 



1 



We then have 



K(m) < 



2(m- LAJ) 



j=-IAI 



IMI 2 



which proves (4.5). 



4.2 Plane wave approximation 

Similar results can be obtained for the plane wave approximation: 

Theorem 4.2 For the space Vf n and the measure v a on the unit disk Q, one has for sufficiently large m 

K(m) > c + cito 2 , (4.12) 
when a = (that is, for the uniform measure), for any c\ < 1/4, where cq depends on c\ and A, and 



K(m) <C\+C 2 



2m 



(4.13) 



when a > 0, for any C% > 1, where C% depends on C2, A and a. 



Proof: Since plane waves are not orthogonal in L 2 (£l, v), the first step is the computation of an orthogonal 
basis of the space spanned by these plane waves. Let us consider 2m + 1 plane waves, with wave vectors 
uniformly distributed on the circle of radius A. For the measures considered here, the Gram matrix of this 
family is a circulant matrix, which is diagonalized in the Fourier basis. An orthogonal family spanning the 
same space is therefore given by functions that are linear combinations of plane waves with the coefficients 
of the discrete Fourier transform: 



1 



(2m + 1) ^ 

j=-m 



e 27rij/(2m+l) ek 



(4.14) 



This formula may be thought as a quadrature for the Bessel integral (2.5 1: the 6™ are thus approximations 
of the Fourier-Bessel functions bj. In order to bound the quantity K(m) = K(V^, v), we compare it to the 



quantity Kiy^v) which behaviour is described by Theorem 4.1 Using (8) and (9) from [5] we have 



itn V^-p(2m+l)t 

°j =2^* '0j+p(2m+l)> 



(4.15) 
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and 



- 1 



* E 

p€Z-{0} 



|ftj+p(2m+l)l 



With m > j > A, we thus have for all < r < 1 and < 9 < 2n, 



bf(r,9) 



- 1 



< 



< 



< 



\bj(r,d)\ 

p6N \ 

\ 2m+l 

7 



+ (p+l)(2m+l 

\P>0 



)M)| +E l & i+(p+l)(2m+l)-2j(r, 
p>0 



-v 



A 

2m+l-2j 



(p+l)(2m+l) 



2m+l 



2A 

3 

1 — A 



where we have used equation (A. 6) of [5] to obtain the second inequality. 
Using the orthogonality of the bj , we have 

W b 7W 2 = E H & j+P(2m+I)l| 2 , 



pez 



and 



|&™|| 2 

- 1 



Nl 2 



E 

pez-{o} 



|&i+p(2m+l) || 2 



IM 2 



When j > A and ( > we bound ||fo J+ i|| 2 according to 



||6, +i || 2 = 2tt / rJ 0+l {Xrfdr 



< 2u r 



21 



Jj(Xr) 2 dr 



< 



21 



2ir / J. } {\rfdr 



21 



where we again have used equation (A. 6) of [5] to obtain the first inequality. We thus have 



WW 2 1 



E H^+p(2m+l)| 
\\bi II 2 



311 \p>0 



\i ; (p+l)(2m+l)l| 2 + l|frj + (p+l)(2m+l)-2jl 
p>0 




We now consider the case a — 0. Using the above comparison results between bj and bj 1 , we can find, 
for any c < 1 a positive integer j c such that for m > j > j c , 



We may thus write 



\b?(l,i 



WbfW 2 



> c 



IIM 2 ' 



(4.19) 



> 



> 



> 



m | 6 „| 2 



E 



3=-m 



I b? II 2 



E 



1^(1,0)1 
||6m||a 



L°°(n) 

2 



v l^(l,0)| a Jj(A)' 



M 2 



INI 2 



The last sum can be bounded from below in a similar way as in the proof of Theorem 4.1 proving (4.12 1. 
We next consider the case a > 0. We then write 



K(m)< J2 



3 = -m 



\h m \\ 2 

imi 2 ' 



(4.20) 



For any C > 1, there is a jc > A such that when j > jc, we have ^\bj(r,9)\ < |6™(r, 6)\ < C|bj(r, so 
that ||6f || L -(n) < C\\b 3 \\ L ^ m = CJj(A), and > |Jj(A)|/C. We then have, for m > j > j c , 



\h m \\ 2 

ii&m 2 



< 



< 



|Lm||2 



l^j\bf\Hx+^J\bf\H6 

D 

1 C 2 Jj(A) 2 
a Jj(X) 2 /C 2 
Ch 
a 



2k 
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and 



K(m) < 



< 



2C 2 (m - [Aj; 



a 



2C 2 (m~ [X\] 



E 

\j\<30 

E 

\j\<jo 



\h m \\ 2 

l°j lli,°°(n) 



IM 2 



which proves 



4.13 



We use here the fact that 116" 



obtained from (|4.15|) and the orthogonality of the bj. 



< 1 which is clear from (4.14), and II 6^ II < \\W 



□ 



5 Numerical tests 

Here, we compare three different reconstructions methods on the unit disc: 

(i) The proposed least-squares method, using Fourier-Bessel approximations. 

(ii) The least-squares method, using the space spanned (2K + l) 2 orthogonal Fourier modes on a square 
enclosing the disc. These modes thus have the form e ak ' x for some fixed < a < ir (here we took 
a = 7r) and k G {— K, . . . , K} 2 and it is assumed that (2K + l) 2 < n. This regularization method is 
one of the most often used in acoustic applications. 

(iii) The Orthogonal Matching Pursuit algorithm used to select a sparse approximation in a dictionary that 
consists of the same Fourier modes, now with (2K + l) 2 > n and K = 40. 

We fix the wave number to be A = 12 and we compare sampling using the measure v a for various proportions 
< a < 1. 



Results for the proposed method (i) are displayed on Figure 5.1 with n = 400 measurements and the 
different values a = 0,0.1,0.5,0.9,1. We plot the error measured in the norm L 2 (il) — L 2 (fl,dx), averaged 
over 10 realizations of the sampling, versus the approximation space dimension 2m + 1. We observe that 
placing more measurement points on the boundary <9f2 is beneficial to the reconstruction: as a approaches 1, 
the lowest error decreases and is attained for larger values of to, confirming our theoretical results. However, 



measuring the solutions only on dfl (a = 1) does not yield good reconstructions. In that case, Theorem 3.1 
only ensures that the reconstruction is accurate on <9f2 and says nothing on the error on the disk itself, since 
we cannot control the L 2 (Vt) norm by the L 2 (n, v{) = L 2 (dil) norm. 



Results with identical parameters for the two other methods (ii) and (iii) are displayed on Figures 5.2 



and 5.3 Reconstruction errors for these methods are larger by a few orders of magnitude than the errors of 
the proposed method (i). Moreover, these two methods do not benefit from sampling on the boundary, as 
the best results are obtained with the uniform sampling. 

Figure [B~4| compares the behavior of these methods "at their best" with varying number of measurements, 
from 100 to 800. The plotted errors are obtained by selecting the value of to for (i), of K for (ii), and of 
the number of iterations of OMP for (iii), as well as the proportion a, that minimize the error for the given 
number n of measurements. Results for (ii) and (iii) are both worse than for the proposed method (i), and 
are always obtained with a = 0. In constrast, the best results of (i) are obtained with a getting closer to 1 
as the number n of measurement grows, for instance a = 0.9 with n = 800 measurements. 

In the practice of using the proposed method, the optimal value of to is unknown to us. A simple way of 
estimating it is by using an independent sample, such as in the generalized cross-validation (GCV) method. 
As an example, for a given number 2to + 1 of Fourier-Bessel functions, we estimate a reconstruction using 
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Figure 5.1: Reconstruction error for method (i) vs. number of Fourier-Bessel functions with n = 400 
measurements 
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Figure 5.2: Reconstruction error for method (ii) vs. number of Fourier modes with n = 400 measurements 




Figure 5.3: Reconstruction error for method (iii) vs. number of OMP iterations, with n = 400 measurements 
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Figure 5.4: Best reconstruction error for method (i) Fourier-Bessel, (ii) Fourier and (iii) OMP 
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Figure 5.5: Best reconstruction error for method (i), using GCV and optimal value of m. 



90% of the samples, and evaluate empirically the mean-square error using the remaining ones. We repeat the 
estimation 10 times using different choices of estimation and reconstruction points within the same sample, 
and select the number m that minimizes the mean square error. The function is then reconstructed using all 
samples. On Figure |5.5| we compare the results obtained by this method with those based on the optimal 
value of m, as the number of measurement n varies. We observe that the performances are comparable up 
to a slight loss by a multiplicative constant. 

6 Conclusion 

In this paper, we compare different ways of sampling solutions of the Hclmholtz equation, using a finite 
number of point measurements. Our main results reveal that good reconstructions can be obtained using 
Fourier-Bessel or plane waves approximations, and that these reconstructions benefit from a denser sampling 
on the boundary. 

These results were obtained in the particular case of a two-dimensional disc. For a more general star- 
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shaped domain in M 2 , the Fourier-Bcssel approximation remains valid. The quantity K(m) does not have an 
explicit expression, yet it can be evaluated numerically after orthogonalization of the Fourier-Bessel family. 
Our first numerical investigations indicate that denser sampling of the functions on the boundary is also 
beneficial in this more general setting. 

In the 3-D case, solutions to the Helmholtz equation in a star-shaped domain can be approximated in 
a similar way, using functions that are products of spherical harmonics and spherical Bcssel functions. We 
expect that a refined analysis of these functions, which have certain common properties with the 2-D Fourier- 
Bessel functions, leads to similar results concerning the convergence of least-squares approximations and the 
role of a denser sampling on the boundary. 
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